CONSTRUCTION OF BAYESIAN DEFORMABLE MODELS 
VIA STOCHASTIC APPROXIMATION ALGORITHM: 
A CONVERGENCE STUDY 
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Abstract. The problem of the definition and the estimation of generative models based on deformable tem- 
plates from raw data is of particular importance for modeling non-aligned data affected by various types of geomet- 
'f0-L ", rical variability. This is especially true in shape modeling in the computer vision community or in probabilistic atlas 

1 building in Computational Anatomy. A first coherent statistical framework modeling the geometrical variability 

^rr' ■ as hidden variables was described by AUassonniere, Amit and Trouve in [2]. The present paper gives a theoretical 

• proof of convergence of effective stochastic approximation expectation strategies to estimate such models and shows 
' the robustness of this approach against noise through numerical experiments in the context of handwritten digit 

(-H modeling. 
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■ 1. Introduction. In the field of image analysis, the statistical analysis and modeling of vaii- 
I able objects from a limited set of examples is still a quite challenging and a largely unsolved 

problem and depends strongly on the use of adequate representations of data. One such rep- 
resentation is the so-called dense deformable template (DDT) framework [4]. Observations are 
defined as deformations, taken from a family of deformations of moderate "dimensionality", of a 
given exemplar or template. Such a representation appears particularly adapted to the emerging 
04 ■ field of Computational Anatomy where one aims at building statistical models of the anatomical 

I variability within a given population [12]. However, research on DDT has been mainly focused 

on the variational point of view, in which DDT is used as an efficient vehicle for a wide range of 
registration algorithms [7]. The problem of template estimation, viewed as a statistical estimation 

■ problem of parameters of generative models of images of deformable objects, has received much 
I less attention. 

' In this paper, we consider the hierarchical Bayesian framework for dense deformable templates 

I developed by AUassonniere, Amit and Trouve in [2] . Each image in a given population is assumed 

' to be generated as a noisy and randomly deformed version of a common template drawn from a 

^ , prior distribution on the set of templates. Individual deformations in their framework are treated 

as hidden variables (or cquivalently random effects in the mixed effects setting), whereas the 
template and the law of the deformations are parameters (or equivalently fixed effects) of interest. 
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■ Parameter estimation for this model could be performed by Maximum A Posteriori (MAP) for 

which existence and consistency (as the number of parameters observed images tends to infinity) 
has been proved (see [2]). This contrasts with earlier work in [11] using a penalized likelihood (PL) 
or the more recent maximum description length approach in [14] for which consistency cannot be 
proved because the deformations are considered as nuisance parameters to be estimated. 

Our contribution in this paper is in defining effective and theoretically proven convergent 
stochastic algorithms for computing (local) maxima of the posterior on the parameters for Bayesian 
deformable template models. First, we specify an adapted stochastic approximation expectation 
minimization algorithm (SAEM algorithm) in this highly demanding framework where the hidden 
variables are non rigid deformation fields living in finite but high dimensional space (typically 
hundreds or more dimensions). In particular, special attention is needed to the sampling of 
the posterior distribution on the deformations. Obviously, MCMC samplers are unavoidable, 
but non adaptive proposal distributions yielding simple symmetric random steps are of limited 
practical interest. The present paper introduces a more sophisticated hybrid Gibbs sampling 
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scheme allowing an acceptable rejection rate during the Estimation step. The overall algorithm 
is cast in the larger class of SAEM-MCMC algorithms introduced in [13]. Second, we extend the 
convergence theory of SAEM-MCMC algorithms developed in [13] to cover the case of unbounded 
random effects arising naturally for deformation fields. The core material of this extension is 
based on the general stability and convergence results for stochastic algorithms with truncation 
on random boundaries given in [6] . The main technical point is that in the presence of unbounded 
random effects and sequential estimation of the covariance matrix of the random effects, the usual 
regularity conditions of the solutions of the Poisson equations for the Markovian dynamic as a 
function of the parameters cannot be verified and have to be relaxed. As a result we provide a new 
general stochastic approximation convergence theorem with a weaker set of assumptions. Third, we 
prove that the conditions for stability and convergence are fulfilled for our general SAEM-MCMC 
estimation algorithm of Bayesian dense deformable templates. Indeed, a well known weakness of 
general stochastic approximation algorithm convergence results is that they rarely provide proofs 
of convergence for the algorithms used in practice since in these implementations the assumptions 
are not satisfied or hard to verify (see [6] ) . Since stochastic approximation algorithms have started 
recently to attract interest for deformable model estimation (see [1] and [16]) our results provide 
the missing theoretical foundations and guidelines for their effective use. As an illustration of the 
potential of such SAEM-MCMC approaches in the context of deformable templates, in particular 
in the presence of noisy data, we present a set of experiments with images of handwritten digits. 

This article is organized as follows. Section 2 briefly reviews the hierarchical Bayesian de- 
formable template model proposed by AUassonnicre, Amit and Trouve in [2]. In Section 3, we 
develop the SAEM-MCMC strategy for the estimation of the parameters. Then in Section 4, we 
state our general convergence result for truncated stochastic approximation algorithm extending 
the Andrieu et al. Theorem of convergence in [6] and state that the designed family of SAEM- 
MCMC algorithms in the previous section satisfy the assumptions. The proof of this last statement 
is postponed to Section 6 after Section 5 concentrates on experiments. In a final Section, we pro- 
vide a short discussion and conclusion. 

2. Observation model. Let us recall the model introduced in [2]. We are given gray level 
images (j/i)i<i<„ observed on a grid of pixels {u„ £ D C G A} which is embedded in a 

continuous domain C M^, (typically D = [—1, 1] x [—1, 1].). Although the images are observed 
only at the pixels («„)„, we are looking for a template image /g : R defined on the plane (the 

extension to images on is straightforward) . Each observation y is assumed to be the discretiza- 
tion on a fixed pixel grid of a deformation of the template plus independent noise. Specifically for 
each observation there exists an unobserved deformation field z : — > such that for m e A 

y(") = loivu - z{vu)) + e(") , 

where e denotes an independent additive noise. 

2.1. Models for template and deformation. Our model takes into account two comple- 
mentary sides: photometry -indexed by p, and geometry -indexed by g. Estimating the template 
and the distribution on deformations directly as a continuous function would be an infinite dimen- 
sional problem. We reduce this problem to a finite dimensional one by restricting the search to a 
parameterized space of functions. The template Iq : — R and the deformation z : R^ ^ R^ 
are assumed to belong to fixed reproducing kernel Hilbert spaces Vp and Vg defined by their re- 
spective kernels Kp and Kg. Moreover, we restrict them to the subset of linear combinations 
of the kernels centered at some fixed control points in the domain D: {vp,j)i<j<k respectively 
i'^g,j)i<3<kg- They are therefore parameterized by the coefficients a G R'^p and f3 G R*"'" x R'^a as 
follows. For all v in D, let 

laiv) = {Kpa){v) =^Kp{v,Vp^.j)a^ , 
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Other forms of smooth parametric representations of the images and of the deformation fields 
could be used without changing the overall results. 

2.2. Parametric model. For clarity, we denote by y* = (?/*,..., y*J and by /3* = . . . , 
the collection of data and their corresponding deformation coefficients. The statistical model of 
the observations we consider is a generative hierarchical one. We assume conditional normal 
distributions for y and (3: 



(2.1) 

/ T 9t T\ \ ^ 9. 

J\A\ 



y ~ «)?=lA/iA|(2ft/a,Cr2ld) I /3,a,cr2 



where denotes the product of distributions of independent variables and zla(u) = Ia(vu — z{vu)), 
for u in A denotes the action of the deformation on the template image. The parameters of 
interest are a which determines the template image, the variance of the additive noise and 
Tg the covariance matrix of the variables (3. We assume that 9 — {a,a'^,Tg) belongs to an open 
parameter space 9: 

&^{e^ {a,<j^,rg) I a e M'^^ j|aj| < i? I, a > 0, Tg G Sym+^ }, 

where ||.|| is the Euclidean norm, Sym^^. is the cone of real positive 2kg x 2kg definite symmetric 
matrices and R an arbitrary positive constant. 

The likelihood of the observed data Qobs can be written as an integral over the unobserved 
deformation variables. Let us denote by the conditional likelihood of the observations given the 
hidden variables and by g„i the likelihood of these missing variables. Then, 



qobs{y\&)= J qc{y\f3,a,<7^)qrnif3\Tg)df3 



where all the densities are determined by the model (2.1). 

2.3. Bayesian model. Even though the parameters are finite dimensional, the maximum- 
likelihood estimator can yield degenerate estimates when the training sample is small. By intro- 
ducing prior distributions on the parameters, estimation with small samples is still possible. The 
regularizing effect of such priors can be seen in the parameter update steps (cf. [2]). We use a gen- 
erative model based on standard conjugate prior distributions for parameters 6 = (a, cr^, Tg) with 
fixed hyper-parameters. Specifically, we assume a normal prior for a, an inverse- Wishart prior on 
and an inverse- Wishart prior on Tg. Furthermore, all priors are assumed to be independent. 
This yields 9 ~ (a, cr^, Fg) ^ Qpara ^ Vp® i^g where 

i'p(da, da^) cc cxp f -i(a - /^p)*(Ep)"^(a - ^p) j fexpf-^j— da'^da, > 3 , 




UgidVg) CC CXpi-{r-\j:g)p/2)^= dV g , flg > 4fcg + 1 . 



(2.2) 

For two matrices A and B, we define {A, B)p = tr{A*B) the Frobenius dot product on the set of 
matrices where tr denotes the trace of the matrix. 

3. Parameter estimation based on stochastic approximation EM. In our Bayesian 
framework, we obtain from [2] the existence of the MAP estimator 

On = argmax(jB(^'|y) , 

see 



4 



S. ALLASSONNIERE, E. KUHN AND A. TROUVE 



where qg denotes the posterior Ukehhood of the parameters given the observations. The depen- 
dence on n refers to the sample size. 

We turn now to the maximization problem of the penalized posterior distribution (7_B(^|y) 
which has no closed form in our ease. Indeed, the probability density function is known up to a 
renormalization constant. That prevents a direct computation of 

In order to solve this problem, we apply an "EM like" algorithm to approximate the MAP 
estimator 6*^. The solution we propose is to base our algorithm on the use of the Stochastic 
Approximation EM (SAEM) . First, we outline certain characteristics of our model, which highlight 
the reasons for the choice of the particular procedure and enable us to simplify its implementation. 

3.1. Model characteristics. An important characteristic of our model is that it belongs to 
the curved exponential family. In other words the complete likelihood q can be written as: 

g(y, (3, 9) = exp + (Si(3), 0(0))] , 

where the sufficient statistic S* is a Borel function on M^, with TV = 2nkg, taking its values in an 
open subset S of M™ and ip, cj) two Borel functions on 8. (Note that S, (/> and ip may depend also 
on y. but since y will stay fixed in what follows, we omit this dependence). 
In our setting, we obtain the following formula: 

logg(y,/3,0) = logq,{y\f3,e) + \ogqrn{f3\e) +logqparai0) , 

where qpara denotes the prior density of the parameters defined in the previous paragraph. 

For any I < j < kp and any u € A, we denote by 

K^{u,j) = Kp{vu - zp{vu),Vpj) 

the matrix which corresponds to the deformation of the kernel Kp through at pixel u and 
evaluated at pixel location 

locpixelu- Then, for some constant C independent of 6, 



logg(y,/3,0) = |-Biog(a2) - ^\\y, K^^aA 



-i logdFgl) - ljr-\j:,)F\ - i(a - Mp)*S;'(a - Mp) 

^2 



c. 



Note that \\yi — K^'aW^ = {yi — K^^aY{yi — K^^a), where K^^a is another way to write the action 
of the deformation zp. on the template denoted previously by z/j.Ia- This form emphasizes 
the dot product between the sufficient statistics and a function of the parameters. It can be 
easily verified that the following matrix-valued functions are the sufficient statistics (up to a 
multiplicative constant) : 



l<i<Tl 
l<i<Tl 



l<i<n 
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For simplicity, wc denote S'(/3) = (5i(/3), S'2(/3), S'3(/3)) for any /3 G and define the sufficient 
statistic space as 

S ^ {(^1, ^2, ^3) I ^1 e ^2 + To^p' e Sym+ , ^3 + agSg G Sym+ ^ } . 

Identifying S2 and S3 with their lower triangular parts, the set S can be viewed as an open set of 
with ns^kp + M|±l) + kg{2kg + 1). 

In [2], the existence of the parameter estimate 0{S) that maximizes the complete log-likelihood 
has been proved. It can easily be shown that a, and V g are explicitly expressed with the above 
sufficient statistics as follows: 

< a{S) = (52+a2(5)(Ep)-i)"'(5i+a2(5)(Ep)-Vp) , (3-1) 
. '^'(^) = ;^A^(«l|y|P + «(5)*52a(5)-2a(5)%+aX) . 
All these formulas also prove the smoothness of 6 on the subset S. 

3.2. SAEM-MCMC algorithm with truncation on random boundaries. In order to 
compute the MAP estimator for our Bayesian model, we use a variant of the EM (Expectation- 
Maximization, [9]) algorithm. This algorithm is quite natural when we have to maximize a likeli- 
hood under a hierarchical model with missing variables. Unfortunately, direct computation is not 
tractable and we have to find a solution to overcome the problematic E step where we have to 
compute an expectation with respect to the posterior distribution on (3 given y. A first attempt 
was proposed in [2] where this conditional distribution is approximated by a Dirac distribution at 
its mode (Fast Approximation with Mode -FAM-EM). The results are very interesting, however, 
the authors point out the lack of convergence of the FAM-EM algorithm when the quality of the 
input images is not good, typically when they are noisy. This is the issue we consider here. We 
propose an algorithm that ensures the convergence of the resulting sequence of estimators toward 
the MAP whatever the quality of the input. 

This solution is a procedure combining the Stochastic Approximation EM (SAEM) with 
Markov Chain Monte Carlo (MCMC) in a more general framework than that proposed by [13], 
which in turn generalized the algorithm introduced by [8]. Indeed, the k*^ iteration of the SAEM- 
MCMC algorithm consists of three steps: 

Step 1 : Simulation step. The missing data, i.e. the deformation parameters /3, are drawn 
using the transition probability of a convergent Markov chain Tig having the posterior 
distribution qpost{-\y-,&) as its stationary distribution: 

/3fc~ne,_,(/3fe_i,.). 

Step 2 : Stochastic approximation step. A stochastic approximation is done on the complete 
log-likclihood using the simulated value of the missing data: 

Qk{0) = Qk-m + Afe_i[logg(y,/3,.,0) - Qu-im , 

where A = (Afc)fe is a decreasing sequence of positive step-sizes. 
Step 3 : Maximization step. The parameters are updated in the M-step: 

Ok = argmax(3fc(6') ■ 

The initial values Qq and ^0 are arbitrarily chosen. 

Remark 1. We cannot use the direct SAEM algorithm. Indeed, this would require sampling 
the hidden variable from the posterior distribution which is known only up to a normalization 
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constant. This sampling is not possible here due to the complexity of the posterior probability 
density function. 

Since our model belongs to the curved exponential family, the stochastic approximation step 
can easily be done on the sufhcient statistics S instead of on the complete log-likelihood. Then 
the maximization step 3 is straightforward, replacing in (3.1) the sufficient statistics with their 
corresponding stochastic approximations. 

The convergence of this algorithm has been proved in [13] in the particular case of missing 
variables living in a compact subset of M^. However, as we set a Gaussian prior on the missing 
variables /3, we cannot assume that their support is compact. In order to provide an algorithm 
whose convergence can be proved in the current framework we have to use a more general setting 
introduced in [6] which involves truncation on random boundaries. The proof is given in Section 
4. This can be formalized as follows. 

Let (/Cg)q>o be a sequence of increasing compact subsets of S such as Ug>o/Cq = S and 
JCq C int(/Cq+i), for all q > 0. Let e = {sk)k>o be a monotone non-increasing sequence of positive 
numbers and K a compact subset of M.^ . We construct a sequence {{(3k, Sk))k>o as described 
in Algorithm 1 as follows. As long as the stochastic approximation does not wander outside 
the current compact set and is not too far from its previous value, we run the SAEM-MCMC 
algorithm. As soon as one of these conditions is not satisfied, we reinitialize the sequences of (3 
and s using a projection (for more details see [6] ), we increase the size of the compact set and 
continue the iterations until convergence. This is detailed in the following steps : 
Initialization step : Initialize Pq and Sg in two fixed compact sets K and ICq respectively. 

Then, for the fc*'* iteration, repeat the following four steps : 

Step 1 : MCMC simulation step . Draw one new element ^ of the non-homogeneous Markov 
Chain with respect to the kernel with the current parameters Ug^, -^ and starting at fB^-i- 



Step 2 : Stochastic approximation step . Compute 

s-s,„i + Ac,_,(^(^)-Sfc-i). (3.2) 

Step 3 : Truncation on random boundaries . If s is outside the current compact set /Ckj..! 

or too far from the previous value Sk , then restart the stochastic approximation in the ini- 
tial compact set, extend the truncation boundary to /C^fc and start again with a bounded 
value of the missing variable. Otherwise, set (/3j.,Sfc) — {(3,s) and keep the truncation 
boundary to JC^k^i ■ 

Step 4 : Maximization step . Update the parameters using (3.1). 

In this algorithm, the MCMC simulation step has to be explained since it involves the choice 
of the transition kernel of the Markov chain. Usually, one uses a Metropolis-Hastings algorithm 
in which a candidate value is sampled from a proposal distribution followed by an accept-reject 
step. However, there arc different possible proposal distributions. The only requirement is that 
all these kernels lead to an ergodic Markov chain whose stationary distribution is our posterior 
distribution. The choice among these possibilities should be based on the specific framework we 
are working in. 

While minimizing the KuUback-Lcibler distance between the stationary distribution /3 —f 
T^eiP) and a tensorial product (3 — > (8i"=iP(/3i) corresponding to independent identically distributed 

n 

missing variables, we get that p is proportional to ^ ^ <lpost{-\yiTO). As n tends to oo and for a 

" 1=1 
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given 9, p converges a.s. towards the prior pdf on the missing variable (7m(-|^)- This suggest to 
use as proposal the prior distribution which involves the current parameters. 

On the other hand, the setting we have in this paper deals with high dimensional missing vari- 
ables. This raises several issues. If we simulate candidates for the hidden variable as a complete 
vector, it appears that most of the candidates are rejected. This is a typical high dimensional 
concentration phenomenon : locally around a current point, the proportion of the space occupied 
by acceptable moves becomes negligible when the space dimension grows. From a more practical 
point of view, even if the proposed candidate is drawn with respect to the current prior distri- 
bution, it creates a deformation that is very different from the current one and too large for the 
corresponding deformed template to fit the observations. This yields very few possible moves 
from the current missing variable value and the algorithm is stuck in a non-optimal location or 
converges very slowly. 

One solution is to update the chain one coordinate at a time conditionally on the others. This 
corresponds to a Gibbs sampler and leads to more relevant candidates which have a higher chance 
to be accepted (cf. [3]). From an image analysis point of view, this put stronger conditions on 
the kind of deformations which are produced when proposing a candidate for each coordinate. 
Knowing the tendency of the movement given by the other coordinates, the candidate will either 
confirm it or not depending if this is a suitable movement. It will thus be accepted with a 
corresponding probability. Even if some coordinates remain unchanged, some others are updated 
which enables the algorithm to visit a larger part of the missing variable support. 



Algorithm 1 Stochastic approximation with truncation on random boundaries 
Set /3q €E K, So G ^o, '>:o = 0, Co = 0, and vq = 0. 
for all A: > 1 do 

compute s = Sk~i + A^^,_^(5(/3) - Sk-i) 

where (3 is sampled from a transition kernel n6/j._j (/3;._]^, .). 
if s e /C„j^_i and |[s - Sfc-i|| < ^Ck-i then 

set (P^, Sk) = (/3, s) and = Kfc-i, i^k = i^fc-i + 1, Cfc = Cfc-i + 1 
else 

set (/3j,, Sk) = s) G K X /Co and Kk = Kk-i + 1, Vk = 0, Cfe = Cfe-i + H^k-i) 
where : N — s- Z is a function such that (f>{k) > —k for any k 
and {f3,s) can be chosen through, different ways (cf . [6]). 
end if 

Ok^ 0{sk). 
end for 



Remark 2. The index k denotes the current active truncation set, the index is the current 
index in the sequences A and e and the index v denotes the number of iterations since the last 
projection. 

3.3. Transition probability of the Markov chain. We now explain how to simulate 
the missing variables thanks to a Markov Chain Monte Carlo algorithm having the posterior 
distribution as its stationary distribution. Due to the inherent high dimensionality N of /3, we 
consider a Gibbs sampler to sequentially scan all coordinates (3-' for I < j < N . 

Denote by (3^^ = {(3^)i^j. We consider here a hybrid Gibbs sampler i.e. each step of the 
Gibbs sampler includes a Metropolis-Hastings step. The proposal law is chosen as qj{-\(3^-' , 6) i.e. 
the conditional law based on the current parameter value 9 derived from the normal distribution 

If 6 is a proposed value at coordinate j, the acceptance rate of the Metropolis-Hastings algo- 
rithm is given by 



q,{b\f3-^,y,9)q,{f3^f3-^,9) 
q,if3^f3-',y,9)qj{b\(3-\9) 



A 1 
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Since 

q, {f3^f3-^ , y, 0) (X q,bs (y 1/3, 0)qj {f3' 1^' , 0) 
the acceptance rate can be simphfied to 



qobs{y\f3,0) 



A 1 



where for any 6 G R and I < j < N , we denote by f3b->j the unique vector which is equal to (3 
everywhere except at coordinate j where it equals b. An illustration of the hybrid Gibbs sampler 
can be found in [17]. The following steps are performed for each coordinate j : 

Step 1 : Proposition . Sample b with respect to the density qj{.\f3^-' ,9). 

Step 2 : Accept-reject . Compute rj{(3\b;fi^\9) and with probability rj{(3^ ,b; (3^^ ,6), up- 
date /3-' to b. 

In Algorithm 2, we summarize the transition step of the Markov chain. 

Algorithm 2 Transition step k ^ k + 1 using a hybrid Gibbs sampler 
Require: f3 ^ f3^.; = 0k 

Gibbs sampler: 

for all j = 1 : N do 

Metropolis-Hastings procedure: 

b^q,{-\(3-\ey, 

Compute r,(/3^ 6; (3~^,9) = [ 'Ztly\m^ ^ ^ 
With probability rj{l3\b\(3^\e), update (3^: (3^ ^ b 
end for 

This yields the transition probability kernel of our Markov chain on /3 : for coordinate j, the 
kernel is 

ne,,(/3,dz) = (®„^,V(dz")) X [<7,(rfz^■|^-^0)r,(/3^dzJ■;/3-^0)+ 



6f3,{dz^) / {l-rj{f3^,b-(3-^,e))q,{b\f3-^,e)db 



(3.3) 

and He = XI^.^v o • • • o Ile.i is therefore the kernel associated with a complete scan. 

4. Convergence analysis. We prove a general theorem on the convergence of stochastic 
approximations for which our algorithm convergence is a special case. 

The hybrid Gibbs sampler used to generate the ergodic Markov chain does not satisfy some 
of the assumptions of the convergence result presented in [6] . We therefore weaken some of their 
conditions, introducing an absorbing set for the stochastic approximation and weakening their 
Holder conditions on some functions of the Markov chain. 

4.1. Stochastic approximation convergence Theorem. Let 5 be a subset of W"-" for 
some integer n^. Let AT be a measurable space. For all s e 5 let Hg : A ^ 5 be a measurable 
function. Let A = {Ak)k be a sequence of positive step-sizes. 

Define the stochastic approximation sequence {sk)k as follows : 

sk = sk-i + Ak-iHs,_,{f3k) with /3fc - ns,_i(/3^,_i, •) , if sfc_i e 5 
Sk = Sc with /3j, = /3c , ii Sk-i S , 

where Sc ^ S, f3^ ^ X and (ns)sg5 is a family of Markov transition probabilities on X. Denote 
by Qa the transition which generates ((/3j., Sfc))^. We consider the natural filtration of the non- 
homogeneous chain {{(3f., Sk))k and denote respectively by and the probability measure 
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and the corresponding expectation generated by this Markov chain starting at (/3, s) and using 
the sequence A. 

If the transition kernel tig of the Markov chain admits a stationary distribution tTs and if for 
any s € 5, Hs is integrable with respect to tTs , then we denote by h the mean field associated with 
our stochastic approximation so that : 

h{s) = J Hsi(3)TTs{f3)df3 . 

The algorithm defined in 4.1 is usually designed to solve the equation h{s) = where h is 
called the mean field function. 

Let {ICq)q>Q be a sequence of increasing compact subsets of S such as U(j>o/Cg = S and 
ICq C int(/Cq+i), Vg > 0. Let e = {ek.)k>o be a monotone non-increasing sequence of positive 
numbers and K a subset of X . 

Let $:Xx5-^Kx/Cobea measurable function and <j> : N Z be a function such that 
4>{k) > —A: for any k. Define the homogeneous Markov chain 

(Zk = (/3fc, Sk,Kk,(k, Vk))k (4.2) 

on Z = X y. S y. 'H^ with the following transition at iteration k : 

• If Vk-i = then draw {(i^^Sk) ^ QA^^_^{^{l3k^ii Sk-i), otherwise draw {f3k,Sk) ~ 

QA(^_^{{f3k-i,Sk-i),-); 

• If \\sk-sk-i\\ < ecfc-i and Sk G /C„^^_l then set Kk = Hk-i, (k = Cfc-i + 1 and i^k = i^k-i + l 
; otherwise set Kk = Kk-i + 1, Cfe = Cfe-i + 4'{^k-i) and Vk = 0. 

Consider the following assumptions, generalized from [6]. Define for any V : X ^ [1, oo] and 
any g : X ^ R"= the norm 

|l.g|ly = sup 



Al'. S is an open subset of M"% h : S is continuous and there exists a continuously 

differentiable function w : 5 ^ [0,oo[ with the following properties. 

(i) There exists an A/o > such that 

£ = {.s e 5, (Vu;(s), /i(s)) = 0} C {s € 5, w{s) < Mq} . 

(ii) There exists a closed convex set 5a C 5 for which s ^ s + pHs{(3) S Sa for any 
p G [0, 1] and (/3, s) € X x Sa {Sa is absorbing) and such that for any Mi g]Mo,oo], 
the set Wmi fl iSa is a compact set of S where Wmi — {s £ 5, w{s) < Mi}. 

(iii) For any s € S\C (V«;(s), h{s)) < 0. 

(iv) The closure of i«(£) has an empty interior. 

A2. For any s G 5, the Markov kernel tig has a single stationary distribution tt^, tTsTIs ~ tTs- In 
addition for all s E S, Hg : X ^ S is measurable and ||7?s(/3)|Ks(d/3) < oo. 

A3'. For any s E S, the Poisson equation g — Ugg = Hs — tTs{Hs) has a solution gg. There exist 
a function V : X ^ [l,oo] such that {/3 E X, V{f3) < oo} ^ 0, constants a g]0, 1], g > 1 
and p > 2 such that for any compact subset /C C 5, 

(i) 

sup \\Hs\\v < oo , (4.3) 
seK 

sup(||5s|k + ||n,g,||y) < oo, (4.4) 



(ii) 



sup ||.s-.s'ir''{|lg, -g.-lly, + |ln,g3-n3,g,,|ly,}<oo. (4.5) 
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(iii) Let fco be an integer. There exist an £ > and a constant C such that for any 
sequence e = {sk)k>o satisfying < < e for all k > k^, for any sequence A ~ 
(Afc)fc>o and for any (3 £ X , 

sup sup E^., [FP«(/3,.)l.(K;)A.(.)>fc] < CVP'^m , (4.6) 
seic fc>o 

where z^(e) = inf{fc > 1, \\sk — Sfc-i|| > and cr(/C) = inf{fc > l,Sfc ^ K.} and the 
expectation is related to the non- homogeneous Markov chain ((/3j., Sk))k>o using the 
step-size sequence {Ak)k>o- 
A4. The sequences A = {Ak)k>o and e = {ek)k>n are non-increasing, positive and satisfy: 

oo oo 

5^ Afe = oo, lim £fc = and J2 {^k + ^fc^fc + i^k^k^T} < where a and p are 

k—O ^ — k—1 

defined in (A3'). 

Theorem 4.1 (General Convergence Result for Truncated Stochastic Approximation). As- 
sume (A1'),(A2), (A3') and (A4). LetKc X be such that sup V{f3) < oo and ICq C Wm,-, n Sa 

(where Mq is defined in (Al')), and let {Zk)k>o be the sequence defined in equation (4.2). Then, 
for all (Sq G K and So £ K-o, we have lim d{sk,C) = OV/s so,o,o,o-o-s, where T/s so. o,o,o 'is the prob- 

ability measure associated with the chain (Zk = {/3j^, Sk, Hk, Ck,i^k))k>a starting at (/3q, so, 0, 0, 0). 

Proof. • The deterministic results obtained by [6] under their assumption (Al) remain true if 
we suppose the existence of an absorbing set as defined in assumption (Al'). Indeed, the proofs in 
[6] can be carried through in the same way restricting the sequences to the absorbing set. Therefore 
we obtain the same properties. The first one (stated in Lemma 2.1 of [6]) gives the contraction 
property of the Lyapunov function w. Then, we have (as in Theorem 2.2 of [6]) the fact that 
a sequence of stochastic approximations stays almost surely in a compact set under some condi- 
tions on the perturbation. Lastly, we establish the convergence of such a stochastic approximation. 

• We then state a relation between the homogeneous and non-homogeneous chains as done in 
Lemma 4.1 of [6]. 

• We now prove an equivalent version of Proposition 5.2 of [6] under our conditions. Indeed 
the upper bound on the fluctuations of the noise sequence stated in this proposition is relaxed in 
our case, involving a different power on the function V. 

Proposition 4.2. Assume (A3'). Let IC he a compact subset of S and let A = {Ak)k and 
£ = {£k)k be two non-increasing sequences of positive numbers such that lim Sk = 0. Then, for p 

k — ^oo 

defined in (A3'), 

1. there exists a constant C such that, for any (/3, s) e X x K, any integer I, any 5 > 



(sup||5i,„(e,A,/C)||><5) <CT-P<| (5]A|) +[Y^Akel] ^ F^'«(/3) 




n>l 



w/iere S'i,„(£, A,/C) = l<T(K;)Ai/(£)>r, ^k{Hs^_i{/3k)-h{sk-i)) andf^^ is the probability 

k=l 

measure generated by the non homogeneous Markov chain {{(3f^,Sk))k started from the 
initial condition (/3,s); 
2. there exists a constant C such that for any {(3, s) G X x /C 



.k=l 
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Proof. The proof of this proposition can proceed as in [6] except for the upper bound on the 
term involving the Holder property (second term in the following). Under (A3'(ii)), this upper 
bound brings into play an exponent pq on the function V. 

Indeed, rewrite 5'i_„(£, A,/C) using the Poisson equation and decompose it into a sum of the 
following five terms : 



fe=l 



.(2) 



k=l 



k=l 
ri-1 

- ^ l::^k^Sk_igsk-Al^k)'^{<j(K)Au(e)=k} ■ 



A5) _ 



(4.7) 

(4.8) 

(4.9) 
(4.10) 
(4.11) 



We evaluate bounds for the first four quantities. Using the Minkowski inequality for p/2 > 1 
and the Burkholder inequality (for T,!^') we have : 



supE^^.s 


sup 

.n>0 


n 


p' 




sup 

.n>0 


n 


P' 


sES 


sup 

.n>0 


n 


P' 


supE^^ 


sup 

.n>0 


n 


P' 



p/2 



<C[Y,^l] sup^E^^,jF''(/3,)l^.(^)^,(,)>fcj] , (4.12) 
<C[f^Aket] sup^E^^^jFf^(/3,)l^,(^)^,(,)>fe}] , (4.13) 



\k=l 



<CA^up^E^^,, [VP{(3,)1 



ses 



{(T(/C)Ay(e)>fc}J : 



(4.14) 



p/2 



a{IC)Av{e)>k}] ■ (4-15) 



\k=l 



ses 



where C is a constant which depends only upon the compact set JC. The higher power pq appears 
because of the Holder condition we assume on the solution of the Poisson equation. 

Since now T',l''^l^o.(^)^^(£)>„} = and noting that V{(3) > 1, V/3 g X, we have V{(3)p < 
VP'^{(3). Applying successively (as in [6]) the Markov inequality, condition (4.6) and the Markov 
property to these upper bounds concludes the proof of the first part of Proposition 4.2. 

Concerning the second part, it follows from the same trick as above for upper-bounding the 
expectation of VP by 

This ends the proof of the proposition. □ 

It is now straightforward to prove the following proposition which corresponds to Proposition 
5.3 in [6]. 

Proposition 4.3. Assume (A3') and (A4)- Then, for any subset K C X such that 

swpV{/3) < oo, any M S (Mo,Afi] and any S > 0, we have lini yl((5, e*"*^', Af , A^'^) = where 
peK fc^oo 

e^'' stands for the sequence e delayed by k switches (ep'^ = £k+i for all I G N) and 

A{S,£,M,A) ^ sup sup (p;f, f sup |15i,fc(e, A, Wm) |1 > s) +P^,(i.(e) < <j{Wm))} ■ 
seJCoPeK L Vfc>i / ' J 
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The convergence of the sequence (sfc)^ follows from the proof of Theorem 5.5 of [6] which states 
the almost sure convergence due to the previous propositions. □ 

Remark 3. We can weaken the condition on p given in (A3'). Indeed, we can assume that 
(A3') holds for any p > as soon as at least condition (4.6) is true also for power 2 ofV. This is 
needed in the proof while giving an upper bound for all the Tn 's using the Jensen inequality instead 
of the Minkowski inequality as in [6]. In this case, the assumption (A4) would have to be satisfied 
for a power max(p, 2) instead of power 2. 

4.2. Convergence Theorem for Dense Deformable Template Model. We now give 
the convergence result of our estimation process which is an application of the previous theorem. 
In this section, we assume that is fixed which reduces 9 to {a, Tg). In fact, due to the implicit 
definition of 9 given in equation (3.1), we were not able to prove the smoothness of the inverse of 
the function s 6{s) which is straightforward for fixed a^. 

We can easily exhibit some of the functions involved in our procedure. Comparing equation 
(3.2) to equation (4.1), we have 

Hs{(3) = S{f3)-s. (4.16) 

Equation (3.1) gives the existence of the hmction s 9{s). We denote by I the observed log- 
likelihood : l{9) = log/ q{y,(3,9)d(3, andlet w(.s) = -loe{s) and /i(.s) = / Hs{l3)qpost{{3\y,0{s))dl3 
for s G 5. 

Theorem 4.4. The sequence of stochastic approximations {sk)k related to the model defined 
in Section 2 and generated by Algorithms 1 and 2 satisfies the assumptions (Al') (ii), (Hi), 
(iv),(A2) and (A3'). 

Proof. The details of the proof are given in appendix (Section 6). □ 

Corollary 1 (Convergence of Dense Deformable Template building via Stochastic Approx- 
imation) . 
Assume 

1. there exist p>l and a e]0, 1[ such that the sequences A — (Afc)i;>o and e = (efc)fc>o are 
non-increasing, positive and satisfy: 

oo QO 

^ Afe = lim Ek ^OandJ2 {^l + ^ksl + (Afce^^)P} < (x; 

k=0 k=l 

2. £ =: {s G 5, (Vi(;(s), /i(s)) = 0} is included in a level set ofw. 

Let K be a compact subset of M.^ and /Co a compact subset of S(M.^). 

Let {sk)k>o o,nd {9k)k>o be the two sequences defined in Algorithms 1 and 2. We denote by 
£' = { 6* e 9{S), ^{9) = 0}, then §{£) = C and 

lim d{9k,C') = P/3o,so,o.o,o - as. 

for all /3q G K and sq G JCq, where P/3q,s(,,o,o,o "is the probability measure associated with the chain 
{Zk = {f3,„Sk,Kk,Ck,i^k))k>o starting at (/3o, sq, 0, 0, 0). 

Proof. We first notice that, as mentioned in [8] (Lemma 2 equation 36), since 9, (p ^'Hd V' are 
smooth functions, it is easy to relate the convergence of the stochastic approximation sequence 
{sk)k to the convergence of the estimated parameter sequence {9k)k- 

Then the proof follows from the general stability result Theorem 4.1 stated in the subsection 
4.1 and from the previous theorem 4.4. □ 

Remark 4. Note that condition (1) is easily checked for Ak = fc^*^ and Sk = k^'^ with 
1/2 < c' < c < 1. However, condition (2) has not been successfully proved yet and should be 
relaxed in future work. 

5. Experiments. To illustrate our stochastic algorithm for the deformable template models, 
we consider handwritten digit images. For each digit class, we learn the template, the correspond- 
ing noise variance and the geometric covariance matrices. (Note that in this experiment the noise 
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variance is no longer fixed and is estimated as the other parameters). We use the US- Postal 
database which contains a training set of around 7000 images. 

Each picture is a (16 x 16) gray level image with intensity in [0, 2] where corresponds to 
the black background. We will also use these sets in the special case of a noisy setting by adding 
independent centered Gaussian noise to each image. 

To be able to compare the results with the previous deterministic algorithm proposed in [2]. 
we use the same samples. In Figure 5.1 below, we show some of the training images. 

OOOS (1000000000600060 

t ! I I M I I I I I I I U I I I I i 

i333 3 5 3 303333^3 33333 
7^1 77717777777777777 

Fig. 5.1. Some images from the training set used for the estimation of the model parameters (inverse video). 




Fig. 5.2. Estimated prototypes of digit 1 (20 images per class) for different hyper-parameters. Left: smoother 
geometry but larger photometric covariance in the spline kernel. Right: more rigid geometry and smaller photo- 
metric covariance. 

M I I ( MHI I I I n M 1 I I 

Fig. 5.3. Synthetic examples corresponding to the two previous estimated templates of digit 1 (inverse video). 
Left : with a fatty shape. Right : with a correct shape thickness. 

A natural choice for the hyper-parameters on a and Vg is ~ Q and we induce the two 
covariance matrices Ep and Eg by the metric of the Hilbert spaces Vp and Vg (defined in Section 
2.1) involving the correlation between the landmarks determined by the kernel. Define the square 
matrices 

Mp{k,k') ^ Kp{vp,j,Vp,j') \fl < k,k' < kp , . 

Mg{k,k') = ifg(%,,%j') VI <k,k'<kg, ^^-'^ 

then Ep = Mp'^ and Eg Mg'^. In our experiments, we have chosen Gaussian kernels for both 
Kp and Kg, where the standard deviations are fixed at ap = 0.12 and ag ~ 0.3. The deformation 
is computed in the [—1,1]^ square with kg = 6 equi-distributed landmarks on this domain. The 
template has been estimated with kp = 15 equi-distributed control points on [—1.5, 1.5]^. 

These two covariance matrices are important hyper-parameters; indeed, it has been shown in 
[2] that changing the geometric covariance has an effect on the sharpness of the template images. 
As for the photometric hyper-parameter, it affects both the template and the geometry in the 
sense that with a large variance, the kernel centered on one landmark spreads out to many of its 
neighbors. This leads to thicker shapes as shown in the left panel of Figure 5.2. As a consequence, 
the template is biased: it is not "centered" in the sense that the mean of the deformations required 
to fit the data is not close to zero. For example for digit "1", the main deformations should be 
contractions or dilations of the template. With a large variance (7^, the template is thicker yielding 
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Fig. 5.4. Estimated prototypes issued from left 10 images per class and right 20 images per class in the 
training set. 

larger contractions and smaller dilations. Since wc have set a Gaussian law on the deformation 
variable P and = —zp, the deformations (Id + zg) and (Id — zp) have the same probability 
to be drawn under the estimated model. As shown on synthetic examples given in Figure 5.3 left 
panel, there are many large dilated shapes. However, these examples were not in the training 
set and are not generated with other hyper-parameters (Figure 5.3 right panel). We have tried 
different relevant values and kept the best with regard to the visual results. We present in the 
following only the results with the adapted variances. 



For the stochastic approximation step-size, we allow a heating period which corresponds to 
the absence of memory for the first iterations. This allows the Markov chain to reach a region of 
interest in the posterior probability density function before exploring this particular region. 

In the experiments presented here, the heating time lasts kh (up to 150) iterations and the 
whole algorithm stops after, at most, 200 iterations depending on the data set (noisy or not). This 
number of iterations corresponds to a point where the convergence seems to have been reached. 
This yields: 

r 1, VI <fc<fc,, 

' ~ I {k\r ' > for d = 0.6 or 1 . 

To optimism the choice of the transition kernel He, we have run the algorithm with different 
kernels and compared the evolution of the simulated hidden variables as well as the results on the 
estimated parameters. Some kernels, as the ones mentioned above do not yield good coverage of 
the infinite support of the unobserved variable. From this point of view the hybrid Gibbs sampler 
we used has better properties and gives nice estimation results which are presented below. 

5.1. Estimated Template. We show here the results of the statistical learning algorithm 
for this model. Figure 5.4 shows two runs of the algorithm for a non-noisy database with 10 and 
20 images per class. Ten images per class are enough to obtain satisfactory template images with 
high contrast. 

Although it was proved in [2] that the KuUback-Leibler divergence between q{-;0) and the 
common density function for observations from a given class converges to its minimal value on 
the family q(-;0), we note that increasing the number of training images does not significantly 
improve the estimated photometric template. This apparently surprising fact can be explained as 
follows : since strong variations in appearance among the images may happened within a given 
class (think about topological changes for instance), the image distribution can not be perfectly 
represented as a distribution around a single template. This distribution is better represented as 
clustered around a major template and minor ones in a multimodal way. When the sample size is 
moderate, with a high probability, the sample contains basically images around the major mode 
and parametric model fits these data quite accurately. When the sample size increases, the minor 
modes start to play a significant role as "outliers" with respect to the major mode in the data, 
resulting in a slightly more blurry template trying to accommodate the different modes. One way 
to overcome this fact is to use some clustering methods as proposed in [2] . To visualize robustness 
with respect to the training set, we ran this algorithm with 20 images per class randomly chosen 
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in the whole database. The different runs are presented in Figure 5.5. The two left images show 
some templates which look like the ones obtained in the left image of Figure 5.4 with the 20 first 
examples of the database. When outliers appear among the 20 randomly chosen training images 
the template may become somewhat more blurry. This is observed for digits '2' and '4' (apparently 
the most variable digits) in the right image of Figure 5.5. For digits where the outliers are less far 
from the other images, the templates are stable. 



mmm u 



smwmai 




Lv] li 





Fig. 5.5. Templates estimated with randomly chosen samples from the whole USPostal database. Each image 
is one run of the algorithm with same initial conditions but different training sets of 20 images per digit each. The 
variability of the results is related to the huge variability inside the USPS database. 
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Fig. 5.6. Evolution of the templates with the algorithm iterations. Top line. Left: Mean gray level images of 
the 20 training samples. Middle: template at the 50th iteration. Right: template at the 100th iteration. Bottom 
line: template at the 150th iteration. The improvement is visible , very fast for some very simple shapes as digit 
1 and longer for very variable ones as digit 2. The higher geometric variability increases the fitting time of the 
algorithm. 

The evolution of the template with the iterations can be viewed in Figure 5.6. The initialization 
of the template is the mean of the gray level images. As the iterations proceed, the templates 
become sharper. In particular, the estimated templates for digits with small geometrical variability 
converge very fast. For digits like '2' or '4', where the geometrical variability is higher, the 
convergence of the coupled parameters (photometry and geometry) is slowed down. 

5.2. Photometric noise variance. The evolution of the noise variance all along the SAEM- 
MCMC iterations is the same as the one observed with the "mode approximation EM" described 
in [2]. During the first iterations, the noise variance balances the inaccuracy of the estimated 
template which is simply the gray-level mean of the training set. As the iterations proceed, the 
template estimates become sharper as does the estimate of the covariance matrix for the geometry. 
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Fig. 5.7. Evolution of the estimated noise variance using 20 images per class along the SAEM-MCMC 
algorithm. This confirms the visual effects seen on the templates : rapid convergence for some really constrained 
digits and slower convergence for very variable ones. 

This yields very small residual noise. Note that here the final noise variance, which is less than 
0.1, for the SAEM-MCMC algorithm for all digits is less than the noise variance , which is between 
0.2 and 0.3, for the mode approximation EM experimented in [2] in the one component run. This 
can be explained by the stochastic nature of the algorithm which enables it to escape from local 
minima provoking early terminations in the deterministic version. 

5.3. Estimated geometric distribution. As mentioned previously we have to fix the value 
of the hyper-parameter ag of the prior on Tg. This quantity plays a significant role in the results. 
Indeed; to satisfy the theoretical conditions we have to choose ag larger than 4A:g -f 1 say 4 x 36-1- 1 in 
our examples. From the geometry update equation, a barycenter between the 'sample' covariance 
and the prior, with the number n of images and ag as coefficients, we find that the prior dominates 
when the training set is small. The covariance matrix stays close to the prior. Thus we need to 
decrease ag and find the best trade-off between the degenerate inverse Wishart and the weight of 
the prior in the covariance estimation. We fix this value with a visual criterion: both the templates 
and the generated sample with the learnt geometry have to be satisfactory. This yields ag = 0.5 
or 0.1. 

As we have observed from Figure 5.9, parameter estimation is robust regardless of whether 
the prior is degenerate or not. In addition, considering the update formulas, even if this law does 
not have a total weight equal to 1 it does not affect parameter estimation. 

In Figure 5.9, we show a sample of some synthetic digits modeled by deformation templates 
drawn with the estimated parameters. Note that the resulting digits in Figure 5.9 look like 
some elements of the training set and seem to explain these data correctly, whereas the prior 
produces some non-relevant local deformations (cf. Figure (5.8)). In particular, for some especially 
geometrically constrained digits such as or 1, the geometry variability reflects their constraints. 
For digits like the 2s, the training set is heterogeneous and shows a large geometrical variability. 
When comparing to the deformations obtained by the mode approximation to EM in [2] , it seems 
that here we obtain a more variable geometry. This might be because with a stochastic algorithm, 
we explore the posterior density and do not only concentrate at its mode. This allows some more 
exotic deformations corresponding to realizations of the missing variable (3 which may belong to 
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Fig. 5.8. Effect of the prior distribution on the deformation : 20 synthetic examples per class generated with 
the estimated template but the prior covariance matrix (inverse video). 
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Fig. 5.9. Effect of the estimated geometric distribution : 40 synthetic examples per class generated with the 
estimated parameters: 20 with the direct deformations and 20 with the symmetric deformations (inverse video). 



the tail of the law. Another reason may be that for such digits, the mode approximation gets stuck 
in a local minimum of the matching energy. Jumping out of this configuration would require a 
large deformation (not allowed by the gradient descent since it would increase the energy again). 
However, such a deformation can be proposed leading to acceptance by the stochastic algorithm. 
Subsequently the deformed template may better fit the observations, leading to acceptance of these 
large deformations. This also leads to a lower value of the residual noise and may also explain the 
low noise variance estimated by the stochastic EM algorithm. 

5.4. Noise effect. As shown in [2], in the presence of noise, the mode approximation al- 
gorithm does not converge towards the MAP estimator. In our setting, the consistency of the 
"SAEM like" algorithm has been proved independently of the training set, and thus noisy images 
can also be treated exactly the same way. These are the results we present here. Figure 5.10 shows 
two training examples per class for noise variance values cr^ = 1 and cr^ = 2. In Figures 5.11 and 
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Fig. 5.10. Two images examples per class of the noisy training set (variance: top: = 1, bottom: = 2). 
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Fig. 5.11. Estimated prototypes in a noisy setting = 1; Left: with the mode approximation algorithm. 
Right: with the SAEM-MCMC coupling procedure. 



5.12, wc show the estimated templates for the noisy training set containing 20 images for both 
methods. Even if the mode approximation algorithm does not diverge, it cannot fit the template 
for digits with a high variability. In contrast, the stochastic EM gives acceptable contrasted tem- 
plates which look like those obtained in Figure 5.4. This becomes more significant as we increase 
the variance of the additive noise we introduce in the training set. 

Concerning the choice of the hyper-parameters, it is not necessary to change all of them. For 
the photometric variance of the spline kernel, a small one could create some non-smooth tem- 
plates and a large kernel would smooth the noise effect. However, we can keep the geometric 
hyper-parameters unchanged. We are presenting here only experiments which seemed to provide 
a reasonable tradeoff between these effects. 



The geometry is also well estimated despite the high level of noise in the training set. Figure 
5.13 shows some synthetic examples, in which parameters are learnt from the training set with 
an additive noise variance of one. The two lines correspond to deformations and their symmetric 
deformation. This sample looks like the synthetic samples learnt on non-noisy images even if some 
examples are not relevant. However, the global behavior has been learnt. 

The algorithm manages to catch the photometry (a contrasted and smoothed template), the 
geometry of the shapes and to "separate" the additive noise. 

The number of iterations needed to reach the convergence point in the noisy setting is about 
twice that of the non-noisy case. The template takes the longest time to converge and the estimate 
of (T^ converges in a few iterations. In particular, the templates obtained in the left panel of Figure 
5.4 with only 10 images per training digit set are obtained with a heating period of 25 iterations 
and 5 more steps with memory. The templates of Figure 5.11, right picture, require 100 to 125 
heating iterations in the 150 global iterations. This is understandable since the algorithm has to 
cope with variations due to the noise and thus needs a longer time to fit the model. 
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Fig. 5.12. Estimated prototypes in a noisy setting = 2: Left: with the mode approximation algorithm. 
Right: with the SAEM-MCMC coupling procedure. 
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Fig. 5.13. Effect of the noise on the geometric parameter estimation : ^0 synthetic examples per class 
generated with the parameters estimated from the noisy training set (additive noise variance ofl, inverse video). 



6. Proof of Theorem 4.4 . Here we demonstrate Theorem 4.4, i.e. the stochastic approxi- 
mation sequence satisfies assumptions (Al') (ii), (iii), (iv), (A2) and (A3'). 

We recah that in this section, the parameter is fixed so that 9 = {a,r). The sufficient 
statistic vector S, the set S as well as the explicit expression of 9{s) have been given in Subsection 
4.2. As noted, ^ is a smooth function of S. 

We will prove that these conditions hold for any p > 1 and a s]0, 1[. 

6.1. Proof of assumption (Al'). 

We recall the functions H, h and w as in [8] defined as follows: 

Hs{(3) = S{f3)-s, 
h{s)^ f Hsif3)qpost{f3\yJis))df3, 

wis) = -ms)). 

As shown in [8], with these functions, we satisfy (Al'(iii)) and (Al'(iv)). 
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Moreover, since the interpolation kernel Kp is bounded, there exist ^4 > and B e Sym^ 
such that for any (3 e M^, we have 

\\Si{fi)\\ < A, < 52(/3) < B and < S^{f3) , 

where, for any symmetric matrices B and B\ we say that B < B' if B' — B is a non-negative 
symmetric matrix. 

We define the set Sa by 

Sa = {SeS \ \\Si\\ <A,0<S2<B and < 5*3 } . 

Since the constraints are obviously convex and closed, we get that Sa is a closed convex subset of 
such that 

SaCSc 

and satisfying 

s + pHsifB) G Sa for any p G [0, 1] any s G Sa and any (3 G M^. 

We now focus on the first two points. As I and 6 are continuous functions, we only need to 
prove that Wm n 5a is a bounded set for a constant M G M.*^_ with: 

Wm = {s G 5, w{s) < M} . 

On Sa, si and S2 are bounded; writing 6{s) — (a(s),r(s)), we deduce from (3.1) and from the 
boundedness of Kp that a{s) is bounded on Sa and \yi — K^'a{s)\ is uniformly bounded on 
Pi G R^'^s and s G 5a. Hence (recall that is fixed here), there exists an 77 > such that 
qc{y\f3, e{s)) > ri for any s G 5a and /3 G M^. Thus, 



W{s) > - log i^J q„^i(3,e{s))df3j+C>-l0g{qparaie{s)))+C>~l0g{qpara\^iris)))+C, 

where C is a constant independent of s G 5'a. Since 

-iog(g,a™|,(r,)) = ^ ((r;\E,)^ + iog|r,|) > ^iog|r,| 



and 



lim log(|rg(s)|) = lim log(|(s3 + agl]g)/(n + %)]) = +00, 

||s||^ + OC, sGOa ||s|l— > + 00,Sf|^5a 



we deduce that 



lim w{s) 



Since w is continuous and Sa is closed, this proves (Al'(ii)). 



6.2. Proof of assumption (A2). 

We prove a classical sufficient condition (DRIl), used in [6] which will imply (A2) under the 
condition that Hs is dominated by V for any s G /C. 

(DRIl) For any s G 5, n^^^^ is ^—irreducible and aperiodic. In addition there exist a function 

V : [1, oo[ and p> 2 such that for any compact subset /C C 5, there exist an integer 

ni and constants 0<A<l,_B>0,K>0,(5>0,a subset C of R^ and a probability 
measure v such that 

supH" F?'(/3) < XV^ifd) + BUf3) , (6.1) 
seic ^ ' 

supn^(^)\/P(/3) < ^VP{(3) V/3 G , (6.2) 

sdlC 

inf nf {j3, A) > Siy{A) V/3 G C, VA G 6(M^) . (6.3) 
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Remark 5. Note that condition (6.3) is equivalent to the existence of a small set C (defined 
below) which only depends on IC. 

Notation 1. Let {ej)i<j<N be the canonical basis o/M^. For any I < j < N, let Eg_j = 
{ /3 G I {f3,ej)g — 0} be the orthogonal space of Span{ej} andpg j be the orthogonal projection 
on Eff j i.e. 



where {f3,/3')g = /^l^g "'^/^i for 9 = {a,Vg) (i.e. the natural dot product associated with the 

covariance matrix Tg) and \\.\\e the corresponding norm. 

We denote for any 1 < j < N and 6 £ Q by T^ej the Markov kernel on M.^ (3.3) associated with 
the Metropolis-Hastings step of the j-th Gibbs sampler step on (3. We have YVg = Ilg jv° ■ ■ 

We first recall the definition of a small set: 

Definition 1. (cf. [15]) A set £ & B{X) is called a small set for the kernel 11 if there exist 
an m > 0, and a non trivial measure Vm on B{X), such that for all (3 £ £ , B £ B{X), 

n"\(3,B)>,,^{B). (6.4) 

When (6.4) holds, we say that £ is Vm-small. 

We now prove the following lemma which give the existence of the small set C in (DRIl): 

Lemma 6.1. Let £ be a compact subset of M.^ and JC a compact subset of S. Then £ is a 
small set o/R^ for n^^^^ for any s fC. 

Proof. First note that there exists an Oc > such that for any G 0, any /3 e and 
any 6 G R, the acceptance rate rj{f3-' , b; f3~'' , 0) is uniformly boimded below by Oc so that for any 
^ l£ j N and any non-negative function /, 

nfl,j/(/3)>ac f fi/3-' +bej)qj{b\f3-',9)db = ac f fipe.ji/3) + ze, /\\ej\\e)9oA^)dz , 

where 50,1 is the density of the standard Gaussian distribution A/'(0, 1). 
By induction, we have 

/ N \ N 

ne/(/3) > / f \pe^NA(^)+'^Zjpe,N,j+i{ej)/\\ej\\e\'[[go,i{zj)dzj, (6.5) 

where pe.q.r — Pe,r ° Pe,r-i o ■ ■ ■ o pg g for any integers q < r and pe,N,N+i = Id. 
Let Ag e £(R^) be the linear mapping on R^ defined by 

AT 

AgZ = ^ZjPg^N,j+liej)/\\ej\\g . 
J = l 

One easily checks that for any 1 < fc < iV, Span{ p0jv.j+i(ej), k < j < N} — Spanjcj \ k < j < 
N} so that Ag is an invertible mapping. By a change of variable, we get 

/ f{pg,N.i{(^) + Aez^)T[9oAzj)dzj^ f{u)gpg ^ ,(p).AeA\{u)du, 

where g^^y, stands for the density of the normal law A/'(/i, S). Since ^ Ag \s, smooth on the set 
of invertible mappings in 0, we deduce that there exist two constants c^; > and Ck. > such 
that CK:Id < AgA^ < Id/cx; and ,9^^ „ (") > C'/c5pfl,«,i(/3),id/cK;(") uniformly for 6 = 9{s) 

with s £ JC. Assuming that f3 £ £, since — > pg,N.i is smooth and £ is compact, we have 
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sup \\pe,N,i{f^)\\ < oo. Therefore, there exist > and c'j^ > such that for any 

pe£fi=e(s), seic 

{u, /3) e X f and any 9 = ^(.s), s e K. 



9pe,N.iiP)^AeAl{u) > C'^goM/c'^i^) ■ (6-6) 

Using (6.5) and (6.6), wc deduce that for any A, for any s G K, and 9 ~ 9{s), 

Uei(3,A)>C'f,a^,y^{A), 

with i/fc equals to the density of the normal law A/'(0, Id/cJ;^). 

This yields the existence of the small set as well as equation (6.3). □ 

This property also implies the </)-irreducibility of the Markov chain (/3j,)fc and its aperiodicity 
(cf. [15] pl21). 

We set V : M.^ [1, +cx3[ as the following function 

V{f3) = l + \\(3\\'- (6.7) 
We, in fact, have the following property : 3 Cic > such that : V/3 G M^, 

sup||i/,(/3)|| <Ck Vi(3). 

This condition is required for the implication of (A2) by (DRIl). 
We now prove condition (6.2). 

Let /C be a compact subset of S and p > 1. For any 1 < j < A^, any s G /C and 9 = 9{s), we 
have 

ne.,VP{(3) < VP{(3) + [ V^{pe.m + ze, /\\e,\\e)g^,,{z)dz . 
Js. 

Since V{l3 + h) < 2{V{l3) + V{h)) for any l3,h & and since there exist two constants cjc > 
and Cjc>0 such that for any /3 G M^, 6* G 9{IC), \\pe,j{l3)\\ < Cjc\\f3\\ and \\ej\\g > 1/cjc, we have 

VP{peA(3) + ze^/\\ej\\e)9oA^^)dz < 2PClVP{f3) jjl + V{c,cze,)YgoAz)dz . 

We deduce that there exists an C'^^ > such that for any (3 G 

sup Ile^,VP{(5)<C'^VP{(3). 

Then, by composition ngVP{l3) < C'f^VP{l3) and (6.2) holds for any p>l. 
Now consider the Drift condition (6.1). 

To prove this inequality, we prove the same inequality for a subsidiary function Vg which 
depends on the parameters 9 and then we deduce the result for V. 
So let us define for any 9 ~ {a, Tg) the function Ve{f3) = 1 + 

Lemma 6.2. Let K be a compact subset of O. For any p > 1, there exist an < px < 1 and 
an Ck > such that for any 9 G K . any (3 G M.^ we have 

UeV,P{(3)<pKVgP{(3) + CK. 
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Proof. The proposal distribution for Tlg^j is given by g(/3 | /3 ,y,9) pgj (/3) + z j^jj^ where 
z ~ M{Q, 1). Then, there exists Ck such that for any f3 G and any measurable set A G S(M^) 

Ug,j{f3,A) ^ {l-ag,i3)lA{l3) + ag,f3 / 1^ ( pg.jil3) + Zj^^) go.i{z)dz , 

where ag_f3 > ac {ac is a lower bound for the acceptance rate). 

Since (pg^j{f3),ej)g = 0, we get Vg (p0,jif3) + zj^'j = Ve{pe^j{f3)) + z^ and 

nejl^J'(/3) - (1 - a8,/3)^/'(/3) + ag^p JjVg{pg^ji(3)) + z'Ygo.Az)dz 

< (1 - ag,i3)VgP{(3) + ag,f, (VgP{pg,,{f3)) + CKVr\poAf3)) ^(1 + zygoAz)dz 

< (1 - ae,/3)y/(/3) + ae,;3F;(pe,,(/3)) + C^.yJ'-^be,,^/?)) • 

We have used in the last inequality the fact that a Gaussian variable has bounded moments of 
any order. Since ag^p > ac and ||p6ij(/3)|je < |j/3|je {pej is an orthonormal projection for the dot 
product (•, •)e), we get that V?7 > 0, 3 Ck^t, such that e and \/ 9 e K 

n,.,FP(/3) < (1 - ac)V^{(3) + {ac + fi)V^{pgAl3)) + Ck,, ■ 

By induction, we show that 

N „ 

He^J'l/S) < ^ - «^)'""^ + (P«.n(/3)) + -^((1 + ^)^+^ - 1) , 

«e{o.i}"i=i ^ 

where pe,u = ((1 - UAr)Id + UNPe,N) o • • • o ((1 — ui)Id + uipe,i). Let = Pe.w ° ■ ■ • ope.i and note 
that pgj is contracting so that 



ngV,^{f3) < bc.^V,P{(3) + {ac + 77)^1/; (pe(/3)) + ^((1 + 77)^+1) , 

for be,, - (E„e{o,i}", nf=i(l - ac)i-"^ K + 77)"^) . 

To end the proof, we need to check that pg is strictly contracting uniformly on K . Indeed, 
lbe(/3)||e = ll/3||e implies that pgj{(3) = f3 for any 1 < j < iV. This yields {/3,ej)g = and thus 
(3 = since {ej)i<j<N is a basis. Using the continuity of the norm of pg in and the compactness 
of K, we deduce that there exists < pK < 1 such that ||pe(/3)|le < Px|l/3||e for any f3 and 6 £ K. 
Changing pk for I > p'j^ > pK we get (1 + pMl3\\lY < p'k{'^ + W/^WeT + C'k for some uniform 
constant C^-. Therefore, 

n.y/(/3) < hc,,V!{^) + p'^l{ac + vfvS{l3) + C'l,,^. 

Since we have inf^>o 6c. r; + p'^k{o-c + v)^ < 1 the result is immediate. □ 
Next, we prove the expected inequality for the function V. 

Lemma 6.3. For any compact set K G Q, any p > 1, there exist < pK < 1, Ck > and 
Too such that Vto > toq , G K , V /3 e 

U'^VP{(3)<pKVP{f3) + CK. 



Proof. Indeed, there exist < ci < C2 such that ciV{f3) < Vg{f3) < C2V{l3) for any {13,9) € 
X K. Then, using the previous lemma, we have U'^'VP{(3) < c7^n^"T//(/3) < cA\PKVe {(3) + 

Ck/{1-Pk)) < {c2/ciy{p^VP{^) + CKl{l- pk))- Choosing to large enough for (c2/ci)fp^ < 1 

gives the result. □ 

This finishes the proof of (6.1) and at the same time of (A2). 
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6.3. Proof of assumption (A3'). 

The geometric ergodicity of the Markov chain, imphed by the Drift condition (6.1), ensures 
the existence of a solution of the Poisson equation (cf [15]): 

%.)(/3) = E(n|(.)^^(/3)-M.'^))- 

A:>0 

We first prove condition (A3'(i)). 

Since Hs{(3) = S{(3) — s with S{j3) at most quadratic in (3, the choice oiV directly ensures (4.3). 

Due to the result presented in [10], there exist upper bounds for the convergence rates and 
the constants involved in the quantification of the geometrical ergodicity of all the chains indexed 
by s €: /C which only depend on m. A, B, S. Therefore, these constants only depend on the fixed 
compact set /C. This yields the uniform ergodicity of the family of Markov chains on /C. So there 
exist constants < 7k; < 1 and Cfc > such that 

\\9e(s)\\v = II E(n|(s)^«(/3) - hismv < ^ C^j^imiv < oo . 

k>Q k>Q 

Thus V.S e /C, gg(^) belongs to Cy ^ {g : R, < oo}. 

Repeating the same calculation as above, it is immediate that Ylg^^^^gg^^^ belongs to Cy too. 
This ends the proof of (A3'(i)). 

We now move to the Holder condition (A3'(ii)). We will use the following lemmas which 
state Lipschitz conditions on the transition kernel and its iterates: 

Lemma 6.4. Let /C fee a compact subset ofS. There exists a constant C/c such that for any 
p > 1 and any function f G Cyp, V(s, s') G IC^ we have : 

||n^(^)/-n^(^,)/|l^,+ ,/. <CK;||/||yp ||.s-,s'|| . 

Proof. For any 1 < j < N and / e Cyp , we have 

ne.,f{(3) = (l-r,(/3,0))/(/3)+ / f{f3,^^)r,{f3\b;f3-\e)q,{b\f3-\9)db, 

where rj(/3, 9) = /g rj{j3-' ,b; f3~-' ,6)qj{b\l3~^ ,6)db is the average acceptance rate. 

Let s and s' be two points in /C and s(e) = (1 — e)s + es' for e S [0, 1] be a linear interpolation 
between s and s' (since S is convex, we can assume that /C is a convex set so that s(e) G JC for 
any e £ [0, 1]). We denote also by 9{e) = ^(s(e)) the associated path in which is a continuously 
differentiable function. To study the difference ||(ne(i)j — ng(o),j)/(/3)|| , introduce Ill jf{f3) = 
(1 - r,(/3,0))/(/3) and Ulj{f3) ^ J^f{f3,_^)r,{f3',b-r^,9)q,{b\f3-^,9)db. We start with the 
difference [[(Ilg^^^^ ^ — Ilg^Q^ j)/(/3)||. First note that under the conditional law qj{b\(3~-' ,9), b ^ 
Uibg^,i/3),l/\\e,\\l) where ' 

beAf3) - e>e,,(/3) = e*/3 - {(3,e,)e/\\e,\\l 
is the j-th coordinate of pg,j{(3). We have 

nL/(/3) - /^/(/3„., +^e,>,(/3^fe;/3-^.)exp 

Since r^ {(3^ , 5; jd'^ , 9) = {(3' ,b;l3-^,9)M where {(3^ , 6; fi-^ , 9) ^ "^^j^iyi^gf ^ is a smooth 
function in 9, we have 

||(n^(,)_^.-n2(o)_^.)/(/3)ll < /V ||/(/3o^,+6e,)| 

Jo JR 



d 



,(/3^fe;/3-^0)exp(- 



{b~beAm\eMM\e 



/27r 



db. 
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However, one easily checks that there exists a constant Cjq such that for any s, s' G /C, e, j and (3 
(with e = 0(e)): 



{b-beAm\e,ro\ \\e,\\o 



— exp 
ae 



<CK;(l + |6-&e.,(/3)|)2exp - 



2tt 

{b-beAm\eM\ \\ej\U 



2tt 



(6.8) 



Since j,\\e,\\e = ^e]±V-\,, fV^' = -T-'fTeT-' and ^r, = ^ (see (3.1)), we 



deduce that there exists another constant Cic such that 



^llejlle <Cic\\s-s\\ 



Similarly, updating the constant Cjc, we have^ 

d 



de 



boAP) 



<C^(1 + !|/3||)P' 



(6.9) 



(6.10) 



Now, concerning the derivative of rj(f3-' , b; f3 \9)^ since 



n 



log(f,(/3^5;/3"^0))^-^ 



'"11 -\m 



1=1 

■th 



with j3^ = f3^ b^jj * corresponding to the i image, only one term of the previous sum is nonzero. 
We deduce from the fact that Kp is bounded and from (3.1) that \-^\og{rj{f3-' ,b; jS'-' ,9))\ < 
^'cl^ckl < Cecils — s'll, so that using the fact that fj{f3-',b\(3~-',6) is uniformly bounded for 
e G e{IC), /3 e M'^ and 6 e M, there exists a new constant C'k: such that 



|^f,(/3^6;/3-^0))| <C^||s-s'|| 



Thus, using (6.8), (6.9) and (6.10), we get for a new constant Cic that 



-r,(/3^6;/3-^0)exp - 



de 



{b-beAm\e 



2 

3\\e 



|ej||e 



< Ceil + m\)h' - + |6- 6e„(/3)|)2exp 



{b-beAm\eX\ \\e,\\e 



Since 11/(^)11 < \\f\\vpVP{p) and V{a + b) = l + \\a + b\\^ < 2{V{a) + V{b)), we have ||/(/3o^, + 
bej)\\ < C\\f\\vp{VP{l3„^j) + VP{bej)) with C = 2'^P-\ Hence, there exists an Ck: such that 
V(s,s') e /C2, V 1 < j < iV, V/3 G andVeG [0,1]: 



ll/(A 



- (r,(/3^6;/3-^0)cxp( - 



(fe-5,,,(/3))2||e,||2\ lle.lle 



< c^ii/ikpi/p(/3o^,)(i + mw - 4 < ccWfWv.v'^m + m\w - 4 



where we have used the fact that a Gaussian variable has finite moments of all order. Since 
(1 + 11/311) < (2F(/3))i/2^ we get (updating C'k) that 



(n^(i),^. - n^(o),,)/(/3)|| < cMv.v^^'^'mw 4 



(6.11) 



^Note that the extra factor (1 + ||/3||) appearing in the RHS of 6.10 compared to the RHS of 6.9 alleviate the 
need to show the usual Lipschitz condition yilg^^j/ — TIg^^,^f\\yj,/ < CicWfWvi \\s — s'\\ with q = p. Weaker 
Lipschitz conditions as conditions A3' (ii) of Theorem 4.1 are needed 
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Now, looking at the first term in (6.3), we deduce easily from the previous study for / = /(/3) 
that 

|!(nj(i),,- - ni(o) J/(/3)|! < CcV {f3)'^'\\s' s\\\\f{f3)\\ < Cdf\\v.VP+'/\(3)\W 4 ; (6.12) 
so that adding (6.11) and (6.12), we get (again updating Cic) that 

IKn^d),, - n0(o),,)/llyp+i/2 < CkW/WvpIW - 4 . (6.13) 

We end the proof, saying that ne(i) - Ugi^o) = Y.j'^i ^e{i),j+i,N ° {'^e{i)j - ^e{Q),j) ° n0(o),ij_i 
where Tlg.q.r = Tlg^r o ne,r-i o • • • o He^q for any integer q < r and any 6* G 8 so that using (6.2) 
and (6.13), the result is straightforward. □ 

Lemma 6.5. Let /C fee a compact subset of S. There exists a constant Cjc such that for all 
p > I and any function f £ Cyp, V(s, s') G K.^ , \fk > 0, we have for 9 = 9{s) and 6' = O(s') that: 

||n^/-n^-,/||v.+v. <Ck||/||v.P-s'|1 . 

Proof. We use the same decomposition of the difference as previously: 

fe-i 

u'gf n^- / = n^(n« - n,0(n^T'-V - ^^^'(/)) ■ 
1=1 

Using Lemma 6.4, the fact that ||ng(/ — 7re(/))||yp < 7^i|/||yp with < 1 (geometric ergodicity) 

and sup sup jjllgy^lly? < cx) wc get: 

j>o eeK 

k-l 

||n^/ - n^,/||y.+i/. <c^Yl " ^e')(^,^-'-V - 7ve'{f))\\v.+u^ 

i=l 

i=l 

and the lemma is proved. □ 

We now prove that h is a, Holder function, adapting linearly Appendix B of [6]. 
Let /3 e and denote by 61 = e{s) and 9' = 9{s'). Write h{s) - h{s') = A{s, s') + B{s, s') + 
C{s, s'), where 

^(,s, s') = {h{s) - U'gHsm + {U'e,HAf3) - his')) , 
i3(s,s') = n^-i/,,(/3)-n^-,i/,(/3), 
Cis,s')^U'g,H,{(3)-U'^,,HAf3). 

Using the geometric ergodicity. Lemma 6.4 and Lemma 6.5, we get that there exists an C > 0, 
independent of k such that: 

\\A{s, s')\\ < C^'' snp \\H,\\vV{l3), 
seic 

\\B{s, s')\\ < C sup \\Hs\\v\\s - s'\\V^/\(5), 
seK 

\\C{s,s')\\<C SUV \\H,\\v\\s-s'\\V{f3). 

This yields 

\\h{s)-h{s')\\<CV'/'mj'' + \\s-s'\\). 
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Hence, setting k — [log \\s — s'||/ log(7)] if \\s — s'\\ < 1 and 1 otherwise, we get the result. 

We can now end the proof of (A3'(ii)): On one hand we have: 

im'^MfB) - his)) - {u'e,HAf3) - his'))\\ < \\u'eHs{f3) - U^gHAM 

+ \\n',HAl3) W^^HAm + \\h{s) h{s')\\ < C\\s s'\\V'/^(3) , 

On the other hand, we have thanks to the geometric ergodicity, 

\m',Hs{(3) - h{s)) - {n'e,HAl3) - h{s'))\\ < Cj'^V'^'iP) . 

Hence for any t > and T > we have 

oc 

||n*g,-(,)(/3) - n^,.gg(,,)(/3)|| < ^ ||(n^i?.(/3) - h{s)) - {n's,HA0) - h{s'j)\\ < 

k=t 



CV^/^ifB) 



T\\s-s'\ 



1-7 



Setting T = [log \\s — s'||/ log(7)] for ||s — s'|| < S < 1 and T = t otherwise, using also the fact that 
for any < a < 1 we have ||s — s'|| log||s — s'\\ = o{\\s — s'||"), we get the result. 

This proves condition (A3'(ii)) for any a < 1. 

We finally focus on the proof of (A3'(iii)). Once again we first prove a specific result for each 
function Ve and obtain after a result for the function V. 

Lemma 6.6. Let K. be a compact subset of S and p > I. There exists Cfz.p > such that for 
any s, s' E JC, for any f3 € M.^ , 

I%)(/3)-^/(,)(/3)I<^^^,pP--^'II%)(/3)- 

Proof. Indeed, there exists C > such that for any 9{s) — {a,Tg) and 9{s') ~ {a',T'g), 
- ^'g\ < C\\s - s'\\. Therefore, there exists an C such that V(s,s') G /C^, [r-^ - {T'g)-^\ < 
C||.s-.s'j| and 

n 

\Veis)i(3) Vgi.,^m < 5]/3*(r, ^ - (r;)-^)A < C\\s - s'\\V{0) . 

i=l 

The result follows from the existence of a constant C such that ^V{(3) < Vg^^j{(3) < CF(/3) for 
any (/3, s) € x /C. □ 

Lemma 6.7. Let tC be a compact subset of S and p> \. There exist e > and C > such 
that for any sequence e = (efc)fe>o such that < e for k large enough, any sequence A = (Afc)fc>o 
and any (3 E M.^ , 

supsupE^.JFP(/3fc)l,(K;)A.(.)>fc] < CVP{(3). 

sGKk>0 

Proof. Let Khe a. compact subset of O such that 9{JC) C K. We note in the sequel, dk = d{sk)- 
We have for fc > 2, using the Markov property and Lemmas 6.2 and 6.6, 

<p{E^jVl_^i(3,_,)K^,c)A.is)>k]+Kti(K^^ +C 

a{K)Aiy{e)>k-l] + ^^'efc-lE^.s [^/^.^ (/3fc-l)l<T(/C)Ai/(£)>fe-l] ) + C . 
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By induction, we show that 



k-l 



C 



(1 _ p(i + C's)) ■ 



Choosing e such that p{l + C'e) < 1 and introducing again < ci < C2 such that ciV{f3) < 
VeifS) < C2V{f3) for any (/3, 0) gR'^ x K end the proof. □ 
This yields (A3'(iii)). 

This concludes the demonstration of Theorem 4.4. 

7. Conclusion and discussion. We have proposed a stochastic algorithm for constructing 
Bayesian non-rigid deformable models in the same context as [2] together with a proof of conver- 
gence toward a critical point of the observed likelihood. To the best of our best knowledge, this is 
the first theoretical result on convergence in the context of deformable template. The algorithm is 
based on a stochastic approximation of the EM algorithm using an MCMC approximation of the 
posterior distribution and truncation on random boundaries. Although our main contribution is 
theoretical, the preliminary experiments presented here on the US-postal database show that the 
stochastic approach can be easily implemented and is robust to noisy situations, yielding better 
results than the previous deterministic schemes. 

Many interesting questions remain open. One may ask what is the convergence rate of such 
stochastic algorithms. A first result has been proved in [8] for the standard SAEM algorithm. 
Under mild conditions, the authors state a central limit theorem for an average sequence of the 
estimated parameters {Ok)k- Concerning the generalization when introducing MCMC, a first step 
has been tackled in [5] . Under some restrictive assumptions the authors can prove a central limit 
theorem for an ergodic adaptive Monte Carlo Markov chain. We truly think that it is possible to 
obtain this kind of convergence rates for the SAEM-MCMC algorithm proposed in this paper. 

Another question refers to the extension of the stochastic scheme to mixture of deformable 
models (defined as the multicomponent model in [2]) where the parameters are the weights of 
the individual components and for each component, the associated template and deformation law. 
This is of particular importance for real data analysis where the restriction to a unique deformable 
model could be too limiting. The design of such mixtures corresponds to some kind of deformation 
invariant clustering approach of the data which is a basic issue in any unsupervised data analysis 
scheme. This extension is, however, not as straightforward as it would appear at first glance: 
due to the high dimensional hidden deformation variables, a naive extension of the Markovian 
dynamics to the component variables will have extremely poor mixing properties leading to an 
impractical algorithm. A less straightforward extension involving multiple MCMC chains is under 
study. 

Another interesting extension is to consider diffeomorphic mappings and not only displacement 
fields for the hidden deformation. This appears to be particularly interesting in the context 
of Computational Anatomy where a one to one correspondence between the template and the 
observation is usually needed and cannot be guaranteed with linear spline interpolation schemes. 
This extension could be done in principle using tangent models based on geodesic shooting in the 
spirit of [18]. 
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